Линейная алгебра. Практическая работа №11

Программные средства

LaTeX. Mathematics

Sage. Linear Algebra

Numpy. Quickstart tutorial

R Functions - Tutorialspoint

Quick-R: User-Defined Functions

Поверхности второго порядка

Общий вид

$F(x,y,z) = 0$

Параметрические

$(t_1, t_2)$ - внутренние координаты поверхности

$x = X(t_1, t_2); y = Y(t_1, t_2); z = Z(t_1, t_2)$

Как функция двух переменных

$z = f(x,z)$

Эллиптический тип

Сфера

$x^2 + y^2 + z^2 = R^2$

Параметрические

$x = R \cos t_1 \sin t_2; y = R \sin t_1 \sin t_2; z = R \cos t_2$

$0 \leq t_1 < 2\pi, 0 \leq t_2 \leq \pi$

Эллипсоид

$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = 1$

Мнимый эллипсоид

$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = -1$

Мнимый конус

$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = 0$

In [1]:
# sage 
x,y,z,R,p,q,r,t1,t2 = var('x,y,z,R,p,q,r,t1,t2')

eq1 = x^2+y^2+z^2-R^2==0
eq2 = x==R*cos(t1)*sin(t2)
eq3 = y==R*sin(t1)*sin(t2)
eq4 = z==R*cos(t2)
eq5 = x^2/p^2+y^2/q^2+z^2/r^2==1
In [2]:
solve([eq1, R==4], R,z)
Out[2]:
[[R == 4, z == -sqrt(-x^2 - y^2 + 16)], [R == 4, z == sqrt(-x^2 - y^2 + 16)]]
In [3]:
solve([eq2, eq3, eq4, R==3], R,x,y,z)
Out[3]:
[[R == 3, x == 3*cos(t1)*sin(t2), y == 3*sin(t1)*sin(t2), z == 3*cos(t2)]]
In [4]:
solve([eq5, p==1,q==2,r==3], p,q,r,z)
Out[4]:
[[p == 1, q == 2, r == 3, z == -3/2*sqrt(-4*x^2 - y^2 + 4)], [p == 1, q == 2, r == 3, z == 3/2*sqrt(-4*x^2 - y^2 + 4)]]
In [5]:
print(sorted(colormaps))
[u'Accent', u'Accent_r', u'Blues', u'Blues_r', u'BrBG', u'BrBG_r', u'BuGn', u'BuGn_r', u'BuPu', u'BuPu_r', u'CMRmap', u'CMRmap_r', u'Dark2', u'Dark2_r', u'GnBu', u'GnBu_r', u'Greens', u'Greens_r', u'Greys', u'Greys_r', u'OrRd', u'OrRd_r', u'Oranges', u'Oranges_r', u'PRGn', u'PRGn_r', u'Paired', u'Paired_r', u'Pastel1', u'Pastel1_r', u'Pastel2', u'Pastel2_r', u'PiYG', u'PiYG_r', u'PuBu', u'PuBuGn', u'PuBuGn_r', u'PuBu_r', u'PuOr', u'PuOr_r', u'PuRd', u'PuRd_r', u'Purples', u'Purples_r', u'RdBu', u'RdBu_r', u'RdGy', u'RdGy_r', u'RdPu', u'RdPu_r', u'RdYlBu', u'RdYlBu_r', u'RdYlGn', u'RdYlGn_r', u'Reds', u'Reds_r', u'Set1', u'Set1_r', u'Set2', u'Set2_r', u'Set3', u'Set3_r', u'Spectral', u'Spectral_r', u'Wistia', u'Wistia_r', u'YlGn', u'YlGnBu', u'YlGnBu_r', u'YlGn_r', u'YlOrBr', u'YlOrBr_r', u'YlOrRd', u'YlOrRd_r', u'afmhot', u'afmhot_r', u'autumn', u'autumn_r', u'binary', u'binary_r', u'bone', u'bone_r', u'brg', u'brg_r', u'bwr', u'bwr_r', u'cool', u'cool_r', u'coolwarm', u'coolwarm_r', u'copper', u'copper_r', u'cubehelix', u'cubehelix_r', u'flag', u'flag_r', u'gist_earth', u'gist_earth_r', u'gist_gray', u'gist_gray_r', u'gist_heat', u'gist_heat_r', u'gist_ncar', u'gist_ncar_r', u'gist_rainbow', u'gist_rainbow_r', u'gist_stern', u'gist_stern_r', u'gist_yarg', u'gist_yarg_r', u'gnuplot', u'gnuplot2', u'gnuplot2_r', u'gnuplot_r', u'gray', u'gray_r', u'hot', u'hot_r', u'hsv', u'hsv_r', u'jet', u'jet_r', u'nipy_spectral', u'nipy_spectral_r', u'ocean', u'ocean_r', u'pink', u'pink_r', u'prism', u'prism_r', u'rainbow', u'rainbow_r', u'seismic', u'seismic_r', u'spring', u'spring_r', u'summer', u'summer_r', u'tab10', u'tab10_r', u'tab20', u'tab20_r', u'tab20b', u'tab20b_r', u'tab20c', u'tab20c_r', u'terrain', u'terrain_r', u'winter', u'winter_r']
In [6]:
g = Graphics()
p1 = implicit_plot3d(lambda x,y,z: x^2+y^2+z^2-4^2, (-4,4), (-4,4), (-4,4), 
                     opacity=0.2, color='#3636ff')
text1 = text3d('x^2+y^2+z^2=4^2', (4,4,4), color='#3636ff', fontsize=12)
p2 = parametric_plot3d((3*cos(t1)*sin(t2),3*sin(t1)*sin(t2),3*cos(t2)), 
                       (t1,0,2*pi), (t2,0,pi), opacity=0.2, color='#ff3636')
text2 = text3d('x=3cos(t1)*sin(t2), y=3sin(t1)*sin(t2), z=3cos(t1)', (-4,-4,4), 
               color='#ff3636', fontsize=12)
p3 = implicit_plot3d(x^2/1^2+y^2/2^2+z^2/3^2-1, (x,-4,4), (y,-4,4), (z,-4,4), 
                     contour=0, color=((z).function(x,y,z),colormaps.bwr_r))
text3 = text3d('x^2/1^2+y^2/2^2+z^2/3^2=1', (-4,4,4), color='#ff36ff', fontsize=12)
g = g+p1+text1+p2+text2+p3+text3
g.show()
In [7]:
# sympy
import sympy
sympy.init_printing(use_unicode=True)
x,y,z,R,p,q,r,t1,t2 = sympy.symbols('x,y,z,R,p,q,r,t1,t2')

eq1 = x^2+y^2+z^2-R^2
eq2 = x-R*cos(t1)*sin(t2)
eq3 = y-R*sin(t1)*sin(t2)
eq4 = z-R*cos(t2)
eq5 = x^2/p^2+y^2/q^2+z^2/r^2-1
In [8]:
sympy.solve([eq1, R-4], R,z)
Out[8]:
$$\left [ \left ( 4, \quad - \sqrt{- x^{2} - y^{2} + 16}\right ), \quad \left ( 4, \quad \sqrt{- x^{2} - y^{2} + 16}\right )\right ]$$
In [9]:
sympy.solve([eq2, eq3, eq4, R-3], R,x,y,z)
Out[9]:
$$\left \{ R : 3, \quad x : 3 \sin{\left (t_{2} \right )} \cos{\left (t_{1} \right )}, \quad y : 3 \sin{\left (t_{1} \right )} \sin{\left (t_{2} \right )}, \quad z : 3 \cos{\left (t_{2} \right )}\right \}$$
In [10]:
sympy.solve([eq5, p-1,q-2,r-3], p,q,r,z)
Out[10]:
$$\left [ \left ( 1, \quad 2, \quad 3, \quad - \frac{3 \sqrt{- 4 x^{2} - y^{2} + 4}}{2}\right ), \quad \left ( 1, \quad 2, \quad 3, \quad \frac{3 \sqrt{- 4 x^{2} - y^{2} + 4}}{2}\right )\right ]$$
In [335]:
# numpy
import numpy
angle = numpy.linspace(0, 2 * numpy.pi, 64)
T1, T2 = numpy.meshgrid(angle, angle)
X = numpy.cos(T1) * numpy.sin(T2)
Y = numpy.sin(T1) * numpy.sin(T2)
Z = numpy.cos(T2)
In [336]:
import pylab as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(4*X,4*Y,4*Z, color='#3636ff', alpha=0.1)
ax.plot_surface(3*X,3*Y,3*Z, color='#ff3636', alpha=0.2)
ax.plot_surface(1*X,2*Y,3*Z, cmap=mpl.cm.bwr_r, linewidth=1, 
                antialiased=True, rstride=2, cstride=2)

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#3636ff', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([1],[1], linestyle="none", c='#ff3636', marker = 'o')
fake2Dline3 = mpl.lines.Line2D([2],[2], linestyle="none", c='#ff36ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2, fake2Dline3], 
          ['$x^2+y^2+z^2=4^2$', 
           '$x=3cos(t1)*sin(t2), y=3sin(t1)*sin(t2), z=3cos(t1)$', 
           '$x^2/1^2+y^2/2^2+z^2/3^2=1$'], loc=9)
ax.set_xlim(-4,4); ax.set_ylim(-4,4); ax.set_zlim(-4,4)
plt.show()
In [115]:
%%r
X1 <- seq(from=-3.4, to=3.4, by=0.75)
Y1 <- X1
sphere1 <- function(R, X, Y){
    func <- function(x, y) {(-x^2-y^2+R^2)^0.5}
    Z <- outer(X, Y, func); Z[is.na(Z)] <- 0
    list(Z, -Z)
}
Z1 <- sphere1(4, X1, Y1)
print(X1); print(Z1[1])







 [1] -3.40 -2.65 -1.90 -1.15 -0.40  0.35  1.10  1.85  2.60  3.35
[[1]]
           [,1]     [,2]      [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]    [,10]
 [1,] 0.0000000 0.000000 0.9110434 1.765644 2.068816 2.077859 1.797220 1.008712 0.000000 0.000000
 [2,] 0.0000000 1.398213 2.3167866 2.766767 2.969428 2.975735 2.787024 2.356905 1.489127 0.000000
 [3,] 0.9110434 2.316787 2.9631065 3.326785 3.497142 3.502499 3.343651 2.994578 2.372762 1.080509
 [4,] 1.7656444 2.766767 3.3267852 3.654449 3.810184 3.815102 3.669809 3.354847 2.813805 1.858763
 [5,] 2.0688161 2.969428 3.4971417 3.810184 3.959798 3.964530 3.824918 3.523847 3.013304 2.148837
 [6,] 2.0778595 2.975735 3.5024991 3.815102 3.964530 3.969257 3.829817 3.529164 3.019520 2.157545
 [7,] 1.7972201 2.787024 3.3436507 3.669809 3.824918 3.829817 3.685105 3.371572 2.833725 1.888783
 [8,] 1.0087121 2.356905 2.9945784 3.354847 3.523847 3.529164 3.371572 3.025723 2.411949 1.164045
 [9,] 0.0000000 1.489127 2.3727621 2.813805 3.013304 3.019520 2.833725 2.411949 1.574802 0.000000
[10,] 0.0000000 0.000000 1.0805091 1.858763 2.148837 2.157545 1.888783 1.164045 0.000000 0.000000

In [114]:
%%r
T1 <- seq(from=0, to=2*pi, length.out=12)
T2 <- seq(from=0, to=pi, length.out=12)
sphere2 <- function(R, t1, t2){list(R*cos(t1)*sin(t2), R*sin(t1)*sin(t2), R*cos(t2))}
sphere2(3, T1, T2)


[[1]]
 [1]  0.000000e+00  7.110255e-01  6.737709e-01 -3.226631e-01 -1.787047e+00 -2.849180e+00
 [7] -2.849180e+00 -1.787047e+00 -3.226631e-01  6.737709e-01  7.110255e-01  3.673940e-16

[[2]]
 [1]  0.000000e+00  4.569484e-01  1.475353e+00  2.244171e+00  2.062362e+00  8.365948e-01
 [7] -8.365948e-01 -2.062362e+00 -2.244171e+00 -1.475353e+00 -4.569484e-01 -8.998559e-32

[[3]]
 [1]  3.0000000  2.8784789  2.5237606  1.9645822  1.2462450  0.4269445 -0.4269445 -1.2462450
 [9] -1.9645822 -2.5237606 -2.8784789 -3.0000000

In [113]:
%%r
X3 <- seq(from=-0.9, to=0.9, by=0.22)
Y3 <- 2*X3
ellipsoid1 <- function(p, q, r, X, Y) {
    func <- function(x, y) {r*(-x^2/p^2-y^2/q^2+1)^0.5}
    Z <- outer(X, Y, func); Z[is.na(Z)] <- 0
    list(Z, -Z)
}
Z3 <- ellipsoid1(1, 2, 3, X3, Y3)
print(X3); print(Z3[1])







[1] -0.90 -0.68 -0.46 -0.24 -0.02  0.20  0.42  0.64  0.86
[[1]]
           [,1]      [,2]      [,3]     [,4]     [,5]     [,6]      [,7]     [,8]      [,9]
 [1,] 0.0000000 0.0000000 0.0000000 1.091604 1.306292 1.161895 0.3498571 0.000000 0.0000000
 [2,] 0.0000000 0.8226786 1.7128923 2.078461 2.198818 2.116223 1.8029975 1.073313 0.0000000
 [3,] 0.0000000 1.7128923 2.2784205 2.564605 2.663081 2.595303 2.3469129 1.846402 0.6627217
 [4,] 1.0916043 2.0784610 2.5646052 2.821914 2.911701 2.849842 2.6256428 2.189795 1.3509996
 [5,] 1.3062925 2.1988179 2.6630809 2.911701 2.998800 2.938775 2.7219111 2.304344 1.5297059
 [6,] 1.1618950 2.1162231 2.5953035 2.849842 2.938775 2.877499 2.6556355 2.225668 1.4084034
 [7,] 0.3498571 1.8029975 2.3469129 2.625643 2.721911 2.655636 2.4134622 1.930285 0.8694826
 [8,] 0.0000000 1.0733126 1.8464019 2.189795 2.304344 2.225668 1.9302850 1.275617 0.0000000
 [9,] 0.0000000 0.0000000 0.6627217 1.351000 1.529706 1.408403 0.8694826 0.000000 0.0000000

Гиперболический тип

Двуполостный гиперболоид

$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = -1$

Однополостный гиперболоид

$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = 1$

Конус

$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = 0$

In [30]:
# sage 
x,y,z,p,q,r = var('x,y,z,p,q,r')

eq1 = x^2/p^2+y^2/q^2-z^2/r^2==-1
eq2 = x^2/p^2+y^2/q^2-z^2/r^2==1
eq3 = x^2/p^2+y^2/q^2-z^2/r^2==0
In [31]:
solve([eq1, p==1,q==2,r==3], p,q,r,z)
Out[31]:
[[p == 1, q == 2, r == 3, z == -3/2*sqrt(4*x^2 + y^2 + 4)], [p == 1, q == 2, r == 3, z == 3/2*sqrt(4*x^2 + y^2 + 4)]]
In [32]:
solve([eq2, p==3,q==4,r==5], p,q,r,z)
Out[32]:
[[p == 3, q == 4, r == 5, z == -5/12*sqrt(16*x^2 + 9*y^2 - 144)], [p == 3, q == 4, r == 5, z == 5/12*sqrt(16*x^2 + 9*y^2 - 144)]]
In [33]:
solve([eq3, p==2,q==3,r==4], p,q,r,z)
Out[33]:
[[p == 2, q == 3, r == 4, z == -2/3*sqrt(9*x^2 + 4*y^2)], [p == 2, q == 3, r == 4, z == 2/3*sqrt(9*x^2 + 4*y^2)]]
In [34]:
g = Graphics()
p1 = implicit_plot3d(x^2/1^2+y^2/2^2-z^2/3^2+1, (x,-7,7), (y,-7,7), (z,-7,7), 
                     contour=0, opacity=0.5, mesh=True, color='#ff36ff')
text1 = text3d('x^2/1^2+y^2/2^2-z^2/3^2=-1', (-7,7,7), color='#ff36ff', fontsize=12)

p2 = implicit_plot3d(x^2/3^2+y^2/4^2-z^2/5^2-1, (x,-7,7), (y,-7,7), (z,-7,7), 
                     contour=0, opacity=0.3, mesh=True, color='#ff3636')
text2 = text3d('x^2/3^2+y^2/4^2-z^2/5^2=1', (7,7,7), color='#ff3636', fontsize=12)
p3 = implicit_plot3d(x^2/2^2+y^2/3^2-z^2/4^2, (x,-7,7), (y,-7,7), (z,-7,7), 
                     contour=0, opacity=0.7,  mesh=True, color='#3636ff')
text3 = text3d('x^2/2^2+y^2/3^2-z^2/4^2=0', (-7,-7,7), color='#3636ff', fontsize=12)
g = g+p1+text1+p2+text2+p3+text3
g.show()
In [71]:
x,y,z,p,q,r = sympy.symbols('x,y,z,p,q,r')

eq1 = x^2/p^2+y^2/q^2-z^2/r^2+1
eq2 = x^2/p^2+y^2/q^2-z^2/r^2-1
eq3 = x^2/p^2+y^2/q^2-z^2/r^2
In [72]:
sympy.solve([eq1, p-1,q-2,r-3], p,q,r,z)
Out[72]:
$$\left [ \left ( 1, \quad 2, \quad 3, \quad - \frac{3 \sqrt{4 x^{2} + y^{2} + 4}}{2}\right ), \quad \left ( 1, \quad 2, \quad 3, \quad \frac{3 \sqrt{4 x^{2} + y^{2} + 4}}{2}\right )\right ]$$
In [73]:
sympy.solve([eq2, p-3,q-4,r-5], p,q,r,z)
Out[73]:
$$\left [ \left ( 3, \quad 4, \quad 5, \quad - \frac{5 \sqrt{16 x^{2} + 9 y^{2} - 144}}{12}\right ), \quad \left ( 3, \quad 4, \quad 5, \quad \frac{5 \sqrt{16 x^{2} + 9 y^{2} - 144}}{12}\right )\right ]$$
In [74]:
sympy.solve([eq3, p-2,q-3,r-4], p,q,r,z)
Out[74]:
$$\left [ \left ( 2, \quad 3, \quad 4, \quad - \frac{2 \sqrt{9 x^{2} + 4 y^{2}}}{3}\right ), \quad \left ( 2, \quad 3, \quad 4, \quad \frac{2 \sqrt{9 x^{2} + 4 y^{2}}}{3}\right )\right ]$$
In [132]:
import warnings
warnings.simplefilter("ignore", category=RuntimeWarning)

X = numpy.linspace(-7, 7, 128)
Y = numpy.linspace(-7, 7, 128)
X, Y = numpy.meshgrid(X, Y) 
def f1(x,y): return 3/2*numpy.sqrt(4*x**2 + y**2 + 4)
def f2(x,y): return 5/12*numpy.sqrt(16*x**2 + 9*y**2 - 144)
def f3(x,y): return 2/3*numpy.sqrt(9*x**2 + 4*y**2)

Z1 = f1(X,Y); Z2 = f2(X,Y); Z3 = f3(X,Y)
Z1[Z1 < -7] = numpy.NaN; Z1[Z1 > 7] = numpy.NaN
Z2[Z2 < -7] = numpy.NaN; Z2[Z2 > 7] = numpy.NaN
Z3[Z3 < -7] = numpy.NaN; Z3[Z3 > 7] = numpy.NaN
In [134]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(X,Y,Z1, color='#ff36ff', alpha=0.7)
ax.plot_wireframe(X,Y,-Z1, color='#ff36ff', alpha=0.7)
ax.plot_wireframe(X,Y,Z2, color='#ff3636', alpha=0.3)
ax.plot_wireframe(X,Y,-Z2, color='#ff3636', alpha=0.3)
ax.plot_wireframe(X,Y,Z3, color='#3636ff', alpha=0.5)
ax.plot_wireframe(X,Y,-Z3, color='#3636ff', alpha=0.5)

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff36ff', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([1],[1], linestyle="none", c='#ff3636', marker = 'o')
fake2Dline3 = mpl.lines.Line2D([2],[2], linestyle="none", c='#3636ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2, fake2Dline3], 
          ['$x^2/1^2+y^2/2^2-z^2/3^2=-1$', 
           '$x^2/3^2+y^2/4^2-z^2/5^2=1$', 
           '$x^2/2^2+y^2/3^2-z^2/4^2=0$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()
In [117]:
%%r
X1 <- seq(from=-0.9, to=0.9, by=0.22)
Y1 <- 2*X1
hyperboloid1 <- function(p, q, r, X, Y) {
    func <- function(x, y) {r*(x^2/p^2+y^2/q^2+1)^0.5}
    Z <- outer(X, Y, func); list(Z, -Z)
}
Z1 <- hyperboloid1(1, 2, 3, X1, Y1)
print(X1); print(Z1[1])






[1] -0.90 -0.68 -0.46 -0.24 -0.02  0.20  0.42  0.64  0.86
[[1]]
          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]
 [1,] 4.855924 4.522345 4.265489 4.099805 4.036533 4.080441 4.228191 4.469497 4.790240
 [2,] 4.522345 4.162115 3.881495 3.698648 3.628388 3.677173 3.840469 4.104632 4.451741
 [3,] 4.265489 3.881495 3.578938 3.379763 3.302726 3.356248 3.534402 3.819791 4.190561
 [4,] 4.099805 3.698648 3.379763 3.168091 3.085774 3.142992 3.332567 3.633841 4.021791
 [5,] 4.036533 3.628388 3.302726 3.085774 3.001200 3.060000 3.254412 3.562303 3.957272
 [6,] 4.080441 3.677173 3.356248 3.142992 3.060000 3.117691 3.308716 3.611980 4.002049
 [7,] 4.228191 3.840469 3.534402 3.332567 3.254412 3.308716 3.489298 3.778095 4.152590
 [8,] 4.469497 4.104632 3.819791 3.633841 3.562303 3.611980 3.778095 4.046332 4.398045
 [9,] 4.790240 4.451741 4.190561 4.021791 3.957272 4.002049 4.152590 4.398045 4.723643

In [127]:
%%r
X2 <- seq(from=-5, to=5, by=1.1)
Y2 <- 4*X2/3
hyperboloid2 <- function(p, q, r, X, Y) {
    func <- function(x, y) {r*(x^2/p^2+y^2/q^2-1)^0.5}
    Z <- outer(X, Y, func); Z[is.na(Z)] <- 0
    list(Z, -Z)
}
Z2 <- hyperboloid2(3, 4, 5, X2, Y2)
print(X2); print(Z2[1])







 [1] -5.0 -3.9 -2.8 -1.7 -0.6  0.5  1.6  2.7  3.8  4.9
[[1]]
           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]     [,10]
 [1,] 10.671874 9.310985 8.137704 7.243771 6.741249 6.718548 7.180220 8.043286 9.195409 10.542243
 [2,]  9.310985 7.713624 6.247222 5.027701 4.272002 4.236088 4.935698 6.123724 7.573712  9.162120
 [3,]  8.137704 6.247222 4.307616 2.192158 0.000000 0.000000 1.972027 4.126473 6.073622  7.966946
 [4,]  7.243771 5.027701 2.192158 0.000000 0.000000 0.000000 0.000000 1.810463 4.810290  7.051399
 [5,]  6.741249 4.272002 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.013865  6.534099
 [6,]  6.718548 4.236088 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.975620  6.510675
 [7,]  7.180220 4.935698 1.972027 0.000000 0.000000 0.000000 0.000000 1.536591 4.714045  6.986097
 [8,]  8.043286 6.123724 4.126473 1.810463 0.000000 0.000000 1.536591 3.937004 5.946521  7.870479
 [9,]  9.195409 7.573712 6.073622 4.810290 4.013865 3.975620 4.714045 5.946521 7.431166  9.044642
[10,] 10.542243 9.162120 7.966946 7.051399 6.534099 6.510675 6.986097 7.870479 9.044642 10.410998

In [128]:
%%r
X3 <- seq(from=-1.9, to=1.9, by=0.4)
Y3 <- 3*X3/2
cone1 <- function(p, q, r, X, Y) {
    func <- function(x, y) {r*(x^2/p^2+y^2/q^2)^0.5}
    Z <- outer(X, Y, func); list(Z, -Z)
}
Z3 <- cone1(2, 3, 4, X3, Y3)
print(X3); print(Z3[1])






 [1] -1.9 -1.5 -1.1 -0.7 -0.3  0.1  0.5  0.9  1.3  1.7
[[1]]
          [,1]     [,2]     [,3]     [,4]      [,5]      [,6]     [,7]     [,8]     [,9]    [,10]
 [1,] 5.374012 4.841487 4.390900 4.049691 3.8470768 3.8052595 3.929377 4.204759 4.604346 5.099020
 [2,] 4.841487 4.242641 3.720215 3.310589 3.0594117 3.0066593 3.162278 3.498571 3.969887 4.534314
 [3,] 4.390900 3.720215 3.111270 2.607681 2.2803509 2.2090722 2.416609 2.842534 3.405877 4.049691
 [4,] 4.049691 3.310589 2.607681 1.979899 1.5231546 1.4142136 1.720465 2.280351 2.952965 3.676955
 [5,] 3.847077 3.059412 2.280351 1.523155 0.8485281 0.6324555 1.166190 1.897367 2.668333 3.452535
 [6,] 3.805260 3.006659 2.209072 1.414214 0.6324555 0.2828427 1.019804 1.811077 2.607681 3.405877
 [7,] 3.929377 3.162278 2.416609 1.720465 1.1661904 1.0198039 1.414214 2.059126 2.785678 3.544009
 [8,] 4.204759 3.498571 2.842534 2.280351 1.8973666 1.8110770 2.059126 2.545584 3.162278 3.847077
 [9,] 4.604346 3.969887 3.405877 2.952965 2.6683328 2.6076810 2.785678 3.162278 3.676955 4.280187
[10,] 5.099020 4.534314 4.049691 3.676955 3.4525353 3.4058773 3.544009 3.847077 4.280187 4.808326

Параболический тип

Эллиптический параболоид

$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 2 * z$

Гиперболический параболоид

$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 2 * z$

Эллиптический цилиндр

$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 1$

Мнимый эллиптический цилиндр

$\frac{x^2}{p^2} + \frac{y^2}{q^2} = -1$

Гиперболический цилиндр

$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 1$

Параболический цилиндр

$y^2 = 2 * p * x$

In [117]:
# sage 
x,y,z,p,q = var('x,y,z,p,q')

eq1 = x^2/p^2+y^2/q^2==2*z
eq2 = x^2/p^2-y^2/q^2==2*z
eq3 = x^2/p^2+y^2/q^2-1==0
eq4 = x^2/p^2-y^2/q^2-1==0
eq5 = y^2==2*p*x
In [118]:
solve([eq1, p==1,q==2], p,q,z)
Out[118]:
[[p == 1, q == 2, z == 1/2*x^2 + 1/8*y^2]]
In [119]:
solve([eq2, p==2,q==3], p,q,z)
Out[119]:
[[p == 2, q == 3, z == 1/8*x^2 - 1/18*y^2]]
In [120]:
solve([eq3, p==3,q==4], p,q,y,z)
Out[120]:
[[p == 3, q == 4, y == 4/3*sqrt(-x^2 + 9), z == r7], [p == 3, q == 4, y == -4/3*sqrt(-x^2 + 9), z == r8]]
In [121]:
solve([eq4, p==2,q==3], p,q,y,z)
Out[121]:
[[p == 2, q == 3, y == 3/2*sqrt(x^2 - 4), z == r9], [p == 2, q == 3, y == -3/2*sqrt(x^2 - 4), z == r10]]
In [122]:
solve([eq5, p==2], p,y,z)
Out[122]:
[[p == 2, y == -2*sqrt(x), z == r11], [p == 2, y == 2*sqrt(x), z == r12]]
In [123]:
g = Graphics()
p1 = implicit_plot3d(x^2/1^2+y^2/2^2-2*z, (x,-6,6), (y,-8,8), (z,-6,6), 
                     contour=0, opacity=0.5, mesh=True, color='#ff3636')
text1 = text3d('x^2/1^2+y^2/2^2=2*z', (7,7,7), color='#ff3636', fontsize=12)

p2 = implicit_plot3d(x^2/2^2-y^2/3^2-2*z, (x,-6,6), (y,-8,8), (z,-6,6), 
                     contour=0, opacity=0.5, mesh=True, color='#3636ff')
text2 = text3d('x^2/2^2-y^2/3^2=2*z', (-7,-7,7), color='#3636ff', fontsize=12)

g = g+p1+text1+p2+text2
g.show()
In [85]:
g = Graphics()
p3 = implicit_plot3d(x^2/3^3+y^2/4^2-1, (x,-6,6), (y,-6,6), (z,-6,6), 
                     contour=0, opacity=0.3,  mesh=True, color='#36ff36')
text3 = text3d('x^2/3^2+y^2/4^2=1', (-7,7,7), color='#36ff36', fontsize=15)
p4 = implicit_plot3d(x^2/2^3-y^2/3^2-1, (x,-6,6), (y,-6,6), (z,-6,6), 
                     contour=0, opacity=0.5,  mesh=True, color='#36ffff')
text4 = text3d('x^2/2^2-y^2/3^2=1', (-7,-7,7), color='#36ffff', fontsize=15)
p5 = implicit_plot3d(y^2-2*2*x, (x,-6,6), (y,-8,8), (z,-6,6), 
                     contour=0, opacity=0.7, mesh=True, color='#ff36ff')
text5 = text3d('y^2=2*2*x', (7,7,7), color='#ff36ff', fontsize=15)
g = g+p3+text3+p4+text4+p5+text5
g.show()
In [86]:
x,y,z,p,q = sympy.symbols('x,y,z,p,q')

eq1 = x^2/p^2+y^2/q^2-2*z
eq2 = x^2/p^2-y^2/q^2-2*z
eq3 = x^2/p^2+y^2/q^2-1
eq4 = x^2/p^2-y^2/q^2-1
eq5 = y^2-2*p*x
In [87]:
sympy.solve([eq1, p-1,q-2], p,q,z)
Out[87]:
$$\left [ \left ( 1, \quad 2, \quad \frac{4 x^{2} + y^{2}}{8}\right )\right ]$$
In [88]:
sympy.solve([eq2, p-2,q-3], p,q,z)
Out[88]:
$$\left [ \left ( 2, \quad 3, \quad \frac{\left(3 x - 2 y\right) \left(3 x + 2 y\right)}{72}\right )\right ]$$
In [89]:
sympy.solve([eq3, p-3,q-4], p,q,y,z)
Out[89]:
$$\left [ \left ( 3, \quad 4, \quad - \frac{4 \sqrt{- \left(x - 3\right) \left(x + 3\right)}}{3}, \quad z\right ), \quad \left ( 3, \quad 4, \quad \frac{4 \sqrt{- \left(x - 3\right) \left(x + 3\right)}}{3}, \quad z\right )\right ]$$
In [90]:
sympy.solve([eq4, p-2,q-3], p,q,y,z)
Out[90]:
$$\left [ \left ( 2, \quad 3, \quad - \frac{3 \sqrt{\left(x - 2\right) \left(x + 2\right)}}{2}, \quad z\right ), \quad \left ( 2, \quad 3, \quad \frac{3 \sqrt{\left(x - 2\right) \left(x + 2\right)}}{2}, \quad z\right )\right ]$$
In [91]:
sympy.solve([eq5, p-2], p,y,z)
Out[91]:
$$\left [ \left ( 2, \quad - 2 \sqrt{x}, \quad z\right ), \quad \left ( 2, \quad 2 \sqrt{x}, \quad z\right )\right ]$$
In [135]:
X = numpy.linspace(-6, 6, 128)
Y = numpy.linspace(-7, 7, 128)
X, Y = numpy.meshgrid(X, Y) 
def f1(x,y): return 1/2*x**2 + 1/8*y**2
def f2(x,y): return 1/8*x**2 - 1/18*y**2

Z1 = f1(X,Y); Z2 = f2(X,Y)
Z1[Z1 < -6] = numpy.NaN; Z1[Z1 > 6] = numpy.NaN
Z2[Z2 < -6] = numpy.NaN; Z2[Z2 > 6] = numpy.NaN
In [136]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(X,Y,Z1, color='#ff3636', alpha=0.7)
ax.plot_wireframe(X,Y,Z2, color='#3636ff', alpha=0.5)

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff3636', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([1],[1], linestyle="none", c='#3636ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2], 
          ['$x^2/1^2+y^2/2^2=2*z$', '$x^2/2^2-y^2/3^2=2*z$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()
In [157]:
X = numpy.linspace(-6, 6, 128)
Z = numpy.linspace(-6, 6, 128)
X, Z = numpy.meshgrid(X, Z) 
def f3(x): return 4/3*(-x**2 + 9)**0.5
def f4(x): return 3/2*(x**2 - 4)**0.5
def f5(x): return -2*(x)**0.5

Y3 = f3(X); Y4 = f4(X); Y5 = f5(X)
Y3[Y3 < -6] = numpy.NaN; Y3[Y3 > 6] = numpy.NaN
Y4[Y4 < -6] = numpy.NaN; Y4[Y4 > 6] = numpy.NaN
Y5[Y5 < -6] = numpy.NaN; Y5[Y5 > 6] = numpy.NaN
In [159]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(X,Y3,Z, color='#36ff36', alpha=0.3)
ax.plot_wireframe(X,-Y3,Z, color='#36ff36', alpha=0.3)
ax.plot_wireframe(X,Y4,Z, color='#36ffff', alpha=0.7)
ax.plot_wireframe(X,-Y4,Z, color='#36ffff', alpha=0.7)
ax.plot_wireframe(X,Y5,Z, color='#ff36ff', alpha=0.5)
ax.plot_wireframe(X,-Y5,Z, color='#ff36ff', alpha=0.5)

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#36ff36', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([0],[0], linestyle="none", c='#36ffff', marker = 'o')
fake2Dline3 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff36ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2, fake2Dline3], 
          ['$x^2/3^2+y^2/4^2=1$', 
           '$x^2/2^2-y^2/3^2=1$',
           '$y^2=2*2*x$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()
In [211]:
%%r
X1 <- seq(from=-6, to=6, by=1)
Y1 <- 2*X1
paraboloid1 <- function(p, q, X, Y) {outer(X, Y, function(x, y) {(x^2/p^2+y^2/q^2)/2})}
Z1 <- paraboloid1(1, 2, X1, Y1)
print(X1); print(Z1)



 [1] -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
 [1,] 36.0 30.5 26.0 22.5 20.0 18.5 18.0 18.5 20.0  22.5  26.0  30.5  36.0
 [2,] 30.5 25.0 20.5 17.0 14.5 13.0 12.5 13.0 14.5  17.0  20.5  25.0  30.5
 [3,] 26.0 20.5 16.0 12.5 10.0  8.5  8.0  8.5 10.0  12.5  16.0  20.5  26.0
 [4,] 22.5 17.0 12.5  9.0  6.5  5.0  4.5  5.0  6.5   9.0  12.5  17.0  22.5
 [5,] 20.0 14.5 10.0  6.5  4.0  2.5  2.0  2.5  4.0   6.5  10.0  14.5  20.0
 [6,] 18.5 13.0  8.5  5.0  2.5  1.0  0.5  1.0  2.5   5.0   8.5  13.0  18.5
 [7,] 18.0 12.5  8.0  4.5  2.0  0.5  0.0  0.5  2.0   4.5   8.0  12.5  18.0
 [8,] 18.5 13.0  8.5  5.0  2.5  1.0  0.5  1.0  2.5   5.0   8.5  13.0  18.5
 [9,] 20.0 14.5 10.0  6.5  4.0  2.5  2.0  2.5  4.0   6.5  10.0  14.5  20.0
[10,] 22.5 17.0 12.5  9.0  6.5  5.0  4.5  5.0  6.5   9.0  12.5  17.0  22.5
[11,] 26.0 20.5 16.0 12.5 10.0  8.5  8.0  8.5 10.0  12.5  16.0  20.5  26.0
[12,] 30.5 25.0 20.5 17.0 14.5 13.0 12.5 13.0 14.5  17.0  20.5  25.0  30.5
[13,] 36.0 30.5 26.0 22.5 20.0 18.5 18.0 18.5 20.0  22.5  26.0  30.5  36.0
In [212]:
%%r
X2 <- seq(from=-6, to=6, by=1)
Y2 <- 3*X2/2
paraboloid2 <- function(p, q, X, Y) {outer(X, Y, function(x, y) {(x^2/p^2-y^2/q^2)/2})}
Z2 <- paraboloid2(2, 3, X2, Y2)
print(X2); print(Z2)



 [1] -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]  [,7]   [,8]   [,9]  [,10]  [,11]  [,12]  [,13]
 [1,]  0.000  1.375  2.500  3.375  4.000  4.375 4.500  4.375  4.000  3.375  2.500  1.375  0.000
 [2,] -1.375  0.000  1.125  2.000  2.625  3.000 3.125  3.000  2.625  2.000  1.125  0.000 -1.375
 [3,] -2.500 -1.125  0.000  0.875  1.500  1.875 2.000  1.875  1.500  0.875  0.000 -1.125 -2.500
 [4,] -3.375 -2.000 -0.875  0.000  0.625  1.000 1.125  1.000  0.625  0.000 -0.875 -2.000 -3.375
 [5,] -4.000 -2.625 -1.500 -0.625  0.000  0.375 0.500  0.375  0.000 -0.625 -1.500 -2.625 -4.000
 [6,] -4.375 -3.000 -1.875 -1.000 -0.375  0.000 0.125  0.000 -0.375 -1.000 -1.875 -3.000 -4.375
 [7,] -4.500 -3.125 -2.000 -1.125 -0.500 -0.125 0.000 -0.125 -0.500 -1.125 -2.000 -3.125 -4.500
 [8,] -4.375 -3.000 -1.875 -1.000 -0.375  0.000 0.125  0.000 -0.375 -1.000 -1.875 -3.000 -4.375
 [9,] -4.000 -2.625 -1.500 -0.625  0.000  0.375 0.500  0.375  0.000 -0.625 -1.500 -2.625 -4.000
[10,] -3.375 -2.000 -0.875  0.000  0.625  1.000 1.125  1.000  0.625  0.000 -0.875 -2.000 -3.375
[11,] -2.500 -1.125  0.000  0.875  1.500  1.875 2.000  1.875  1.500  0.875  0.000 -1.125 -2.500
[12,] -1.375  0.000  1.125  2.000  2.625  3.000 3.125  3.000  2.625  2.000  1.125  0.000 -1.375
[13,]  0.000  1.375  2.500  3.375  4.000  4.375 4.500  4.375  4.000  3.375  2.500  1.375  0.000
In [313]:
%%r
X3 <- seq(from=-5, to=5, by=0.1)
Z3 <- X3
cylinder1 <- function(p, q, X, Z) {
    for (x in X) {
        for (z in Z) {
            y <- q*(1-x^2/p^2)^0.5
            v <- c(x,y,z,x,-y,z)
            if (sum(is.na(v)) == 0) {
                write.table(matrix(v, nrow=2, ncol=3, byrow=T), 
                            file ="cylinder1.csv", append = TRUE, quote = FALSE,
                            col.names = FALSE, row.names = FALSE)
            }
        }
    }
}

cylinder1(3, 4, X3, Z3)















In [314]:
%%r
Y4 <- seq(from=-5, to=5, by=0.1)
Z4 <- X4
cylinder2 <- function(p, q, X, Z) {
    for (x in X) {
        for (z in Z) {
            y <- q*(x^2/p^2-1)^0.5
            v <- c(x,y,z,x,-y,z)
            if (sum(is.na(v)) == 0) {
                write.table(matrix(v, nrow=2, ncol=3, byrow=T), 
                            file ="cylinder2.csv", append = TRUE, quote = FALSE,
                            col.names = FALSE, row.names = FALSE)
            }
        }
    }
}

cylinder2(2, 3, X4, Z4)















In [315]:
%%r
Y5 <- seq(from=-5, to=5, by=0.1)
Z5 <- Y5
cylinder3 <- function(p, Y, Z) {
    for (y in Y) {
        for (z in Z) {
            x <- y^2/2/p
            v <- c(x,y,z)
            write.table(matrix(v, nrow=1, ncol=3, byrow=T), 
                        file ="cylinder3.csv", append = TRUE, quote = FALSE,
                        col.names = FALSE, row.names = FALSE)
        }
    }
}

cylinder3(2, Y5, Z5)













In [328]:
X1, Y1, Z1 = numpy.loadtxt('cylinder1.csv', delimiter=' ', unpack=True)
X2, Y2, Z2 = numpy.loadtxt('cylinder2.csv', delimiter=' ', unpack=True)
X3, Y3, Z3 = numpy.loadtxt('cylinder3.csv', delimiter=' ', unpack=True)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X1, Y1, Z1, s=1, c='#36ff36', marker = '1')
ax.scatter(X2, Y2, Z2, s=30, c='#36ffff', marker = '3')
ax.scatter(X3, Y3, Z3, s=3, c='#ff36ff', marker = '1')
fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#36ff36', marker = '1')
fake2Dline2 = mpl.lines.Line2D([0],[0], linestyle="none", c='#36ffff', marker = '3')
fake2Dline3 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff36ff', marker = '1')
ax.legend([fake2Dline1, fake2Dline2, fake2Dline3], 
          ['$x^2/3^2+y^2/4^2=1$', 
           '$x^2/2^2-y^2/3^2=1$',
           '$y^2=2*2*x$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()

Плоскости

Пара мнимых пересекающихся плоскостей

$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 0$

Пара пересекающихся плоскостей

$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 0$

Пара параллельных плоскостей

$y^2 = p^2$

Пара мнимых параллельных плоскостей

$y^2 = -p^2$

Пара совпадающих плоскостей

$y^2 = 0$

In [196]:
# sage 
x,y,z,p,q = var('x,y,z,p,q')

eq1 = x^2/p^2-y^2/q^2==0
eq2 = y^2==p^2
eq3 = y^2==0
In [197]:
g = Graphics()
p1 = implicit_plot3d(x^2/1^2-y^2/2^2, (x,-3,3), (y,-3,3), (z,-3,3), 
                     contour=0, color=(sin(x+z)*cos(y+z),colormaps.bwr_r))
text1 = text3d('x^2/1^2=y^2/2^2', (4,4,4), color='#ff3636', fontsize=12)

p2 = implicit_plot3d(y^2-1^2, (x,-3,3), (y,-3,3), (z,-3,3), 
                     contour=0, opacity=0.7, color='#ff36ff')
text2 = text3d('y^2=1^2', (-4,-4,4), color='#ff36ff', fontsize=12)

g = g+p1+text1+p2+text2
g.show()
In [117]:
solve([eq1, p==1,q==2], p,q,y,z)
Out[117]:
[[p == 1, q == 2, y == -2*x, z == r29], [p == 1, q == 2, y == 2*x, z == r30]]
In [118]:
solve([eq2, p==1], p,y,x,z)
Out[118]:
[[p == 1, y == -1, x == r31, z == r32], [p == 1, y == 1, x == r33, z == r34]]
In [119]:
# sympy
x,y,z,p,q = sympy.symbols('x,y,z,p,q')

eq1 = x^2/p^2-y^2/q^2
eq2 = y^2-p^2
eq3 = y^2
In [122]:
sympy.solve([eq1, p-1,q-2], p,q,y,z)
Out[122]:
$$\left [ \left ( 1, \quad 2, \quad - 2 x, \quad z\right ), \quad \left ( 1, \quad 2, \quad 2 x, \quad z\right )\right ]$$
In [123]:
sympy.solve([eq2, p-1], p,y,x,z)
Out[123]:
$$\left [ \left ( 1, \quad -1, \quad x, \quad z\right ), \quad \left ( 1, \quad 1, \quad x, \quad z\right )\right ]$$
In [282]:
X = numpy.linspace(-6, 6, 128)
Z = numpy.linspace(-6, 6, 128)
X, Z = numpy.meshgrid(X, Z) 
def f1(x): return 2*x
def f2(x): return 1
Y1 = f1(X); Y2 = f2(X)
Y1[Y1 < -6] = numpy.NaN; Y1[Y1 > 6] = numpy.NaN
In [283]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(X,Y1,Z, cmap=mpl.cm.bwr_r, linewidth=1, 
                antialiased=True, rstride=2, cstride=2)
ax.plot_surface(X,-Y1,Z, cmap=mpl.cm.bwr_r, linewidth=1, 
                antialiased=True, rstride=2, cstride=2)
ax.plot_wireframe(X,Y2,Z, color='#ff36ff', alpha=0.5)
ax.plot_wireframe(X,-Y2,Z, color='#ff36ff', alpha=0.5)

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#3636ff', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([1],[1], linestyle="none", c='#ff36ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2], 
          ['$x^2/1^2=y^2/2^2$', '$y^2=1^2$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()
In [329]:
%%r
X <- seq(from=-6, to=6, by=0.15)
Z <- X
planes1 <- function(p, q, X, Z) {
    for (x in X) {
        for (z in Z) {
            y <- x*q/p
            if (min(X) <= y && y <= max(X)) {
                v <- c(x,y,z,x,-y,z)
                v <- replace(v, is.na(v), 0)
                write.table(matrix(v, nrow=2, ncol=3, byrow=T), 
                            file ="planes1.csv", append = TRUE, quote = FALSE,
                            col.names = FALSE, row.names = FALSE)
            }
        }
    }
}

planes2 <- function(p, X, Z) {
    for (x in X) {
        for (z in Z) {
            v <- c(x,p,z,x,-p,z)
            write.table(matrix(v, nrow=2, ncol=3, byrow=T), 
                        file ="planes2.csv", append = TRUE, quote = FALSE,
                        col.names = FALSE, row.names = FALSE)
        }
    }
}

planes1(1, 2, X, Z); planes2(1, X, Z)


























In [334]:
X1, Y1, Z1 = numpy.loadtxt('planes1.csv', delimiter=' ', unpack=True)
X2, Y2, Z2 = numpy.loadtxt('planes2.csv', delimiter=' ', unpack=True)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X1, Y1, Z1, s=8, c=Z1, cmap=mpl.cm.bwr_r, marker='1')
ax.scatter(X2, Y2, Z2, s=2, c='#ff36ff', marker='1')

fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#3636ff', marker='1')
fake2Dline2 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff36ff', marker='1')
ax.legend([fake2Dline1, fake2Dline2], 
          ['$x^2/1^2=y^2/2^2$', 
           '$y^2=1^2$'], loc=9)

ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()